Survo C codes for MAT #QFIND, MAT #QRFIND, and CTEST:


/* #qfind.c 3.5.2017/SM (17.5.2017)
*/

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <malloc.h>
#include <math.h>
#include <survo.h>
#include <survoext.h>
#include "ext_mat.h"

extern char *argv1;
static int i,n,h,j,jj;
static int n21,i0;
static int row,col;
static int *Q,*perm;

static int *P,*P2;
extern char **spb;
int np;
static char name[64];

// MAT #QFIND(n) or MAT #QFIND(n,TEST)

FILE *tmp;

op__qfind()
        {

        i=external_mat_init(1); if (i<0) return(1);
        if (g<3)
            {
            init_remarks();
            rem_pr("MAT #QFIND(n,TEST)                            ");
            rem_pr("http://www.survo.fi/papers/Roots2013.pdf ");
            rem_pr("creates nxn Qn table (2017)              ");
            rem_pr("Optional TEST checks the properties of Qn");
            wait_remarks(2);
            return(1);
            }

        n=1+atoi(word[2]);

        np=0;
        if (g>3) np=n;
        if (i>0) { P=malloc(np*sizeof(int)); P2=malloc(np*sizeof(int)); }

        sprintf(sbuf,"Q%d.TXT",n-1);
        tmp=fopen(sbuf,"wt");

        n21=2*n+1;
        Q=malloc(2*n21*(n+1)*sizeof(int));
        perm=malloc(2*(n+1)*sizeof(int));
        for (i=0; i<n21*(n+1); ++i) Q[i]=0;
        for (i=0; i<n; ++i) Q[n21+n+1+i]=i+1;

        row=0; col=2;
        while (row<n-1)
            {
            ++row;
            for (i=0; i<row; ++i) perm[i]=Q[n21+n+row+i*n21];
            if (np)
                {
                for (i=0; i<row; ++i) P[i]=perm[i];
                test(row);
                }
            for (i=0; i<row; ++i) fprintf(tmp,"%d\n",perm[i]);

            set_row(row);
            h=col*n21+n+row;

            j=h-2*col+1; jj=h;
            while (jj-j>2)
                {
                Q[j++]=-Q[jj--];
                }

            j=h+1;
            while (j<(col+1)*n21)
               {
               Q[j]=2*Q[j-col]-Q[j-2*col];
               ++j;
               }
            ++col;
            }

        fclose(tmp);

        external_mat_end(argv1);
        return(1);
        }

set_row(k)
int k;
    {
    int j,h,i;

    j=1;
    h=n21+n+k;
    while (j<=n)
        {
        for (i=0; i<k; ++i) { Q[h]=perm[i]; j++; if (j>n) break; h+=n21; }
        if (j>n) break;
        for (i=0; i<k; ++i) { Q[h]=perm[k-1-i]; j++; if (j>n) break; h+=n21; }
        Q[h]=0; j++; h+=n21; if (j>n) break;
        }
    return(1);
    }

// None of error messages in this function should appear!
int test(n)
int n;
    {
    int i,j;

    if (P[0]!=n)
        { printf("\nQFIND: The first element is not the greatest one %d !", n);
          getch(); exit(0);
        }
    for (i=0; i<n; ++i)
        if (P[i]==0) P[i]=i+1; // Replace 0 by column index
    for (i=0; i<n; ++i) P2[i]=0;
    for (i=0; i<n; ++i)
        {
        j=P[i];
        if (P2[j-1]!=0)
           { printf("\nQFIND: %d appears at least twice!",j);
             getch(); exit(0);
           }
        P2[j-1]=j;
        }

/***************
  Testing: All permutations of order 2 ?
  row=10
  j       0  1  2  3  4  5  6  7  8   cycles (1,9)(2,5)(3)(4,7)(6,8)
          1  2  3  4  5  6  7  8  9
  P[j]    9  5  3  7  2  8  4  6  1
**************/

    for (j=0; j<n; ++j)
        {
        if (P[P[j]-1]!=j+1)
           {
     printf("\nQFIND: Not permutation of order 2: %d %d",P[P[j]-1],j+1);
           getch(); return(-1);
           }
        }
    return(1);
    }


/* #qrfind.c 10.7.2013/SM (24.4.2014) (11.4.2017) (5.5.2017) (21.6.2017)
*/

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <malloc.h>
#include <math.h>
#include <survo.h>
#include <survoext.h>
#include "ext_mat.h"

extern char *argv1;
static int i,n,m,k,j,q,i1;
static double *D,*H;
static double *X;
static double *C;
static int *Q;
static int am;

static FILE *tmp;

static char name[64];
static char *lab;

// MAT #QRFIND(n,Q,matrix)

extern double pi;
extern double sgn();
extern double cot();

op__qrfind()
        {
        char expr1[2*LLENGTH];
        double a,b,gg,gap;
        int odd=0;
        int c_q,c_factor,c_index,c_amod;
        int f;

        i=external_mat_init(1); if (i<0) return(1);
        if (g<4)
            {
            init_remarks();
            rem_pr("MAT #QRFIND(n,Q,A)                                     ");
            rem_pr("http://www.survo.fi/papers/Roots2013.pdf             ");
            rem_pr("Numerical roots directly as cot expressions! (2017)  ");
            rem_pr("                                                     ");
            rem_pr("sign, q, factor, index, amod  values saved as A.MAT  ");
            rem_pr("q values obtained from Q file generated by           ");
            rem_pr("MAT #QFIND.                                          ");
            rem_pr("Matrix of coefficients saved as Cn.MAT               ");
            wait_remarks(2);
            return(1);
            }

        n=atoi(word[2]);
        odd=1; if (n%2==0) odd=0;
        m=n/2;

        i=spec_init(r1+r-1); if (i<0) return(-1);

        Q=(int *)calloc(m*m,sizeof(int));

        tmp=fopen(word[3],"rt");

        for (i=0; i<m; ++i)
            {
            for (j=0; j<=i; ++j)
                 {
                 fscanf(tmp,"%d\n",&Q[j+m*i]);
                 if  (np) P[j]=Q[j+m*i]; // 6.5.2017  test permutation
                 }
            }

        mat_alloc(&X,m,1); // Numerical roots
        for (i=0; i<m; ++i) X[i]=cot((2*i+1)*3.141592653589793/(2*n));
        // X corresponds to R/n.

        odd=n-2*m; k=2-odd;

        c_q=m;
        c_factor=2*m;
        c_index=3*m;
        c_amod=4*m;

        i=mat_alloc_lab(&T,m,5,&rlabT,&clabT);
        strcpy(clabT,"sign    q       factor  index   amod");
        numlab(rlabT,m,8);
        for (i=0; i<5*m; ++i) T[i]=0.0;
        D=malloc(m*sizeof(double));
        H=malloc(m*sizeof(double));

        for (i=1; i<=m; ++i)
            {
            D[i-1]=2*sin((double)(m+1-i)*pi/(double)n);
            T[c_amod+i-1]=(double)amod(n,2*i-1); // 24.4.2014
            }
        if (odd==0) D[0]=1;
        // D corresponds to E/n.

        i=1;
        while (i<=m)
            {
            if (i==1)
                 {
                 if (n%2==0) T[i-1]=1; else T[i-1]=-1;
                 T[c_q+i-1]=1;
                 ++i; iprint(i); continue;
                 }
            if (T[c_factor+i-1]>0) { ++i; q=1; continue; }
            am=amod(n,2*i-1);

            if (am==0) q=0; else q=Q[(i-2)*m+am-1];
            if (q>0)
                {
                if (n%2==0)
                   {
                   if (i%2==0) T[i-1]=-1; else T[i-1]=1;
                   }
                else
                   {
                   if (q%2==i%2) T[i-1]=-1; else T[i-1]=1;
                   }
                T[c_q+i-1]=q;
                ++i;
                continue;
                }
            else
                {
                // If no valid q is found,          Q
                // X[i-1] must be a 'factor' root.
                a=X[i-1];
                f=(int)(pi/(2*atan(1.0/a))+0.5);
                if (n%f!=0)
                    {
                    // This error message should never appear!
                    sprintf(sbuf,"\n%d is not a factor of %d!",
                                  f,n);
                    sur_print(sbuf); getch();
                    return(1);
                    }
                gap=n/f; // Recording other roots related to f
                i1=i; q=1;
                while (i1<=m)
                    {
                    T[c_factor+i1-1]=f;
                    T[c_index+i1-1]=q;
                    i1+=(int)gap; ++q;
                    }
                ++i; iprint(i); q=1; continue;
                }
            } // i

        mT=m;
        nT=5;
        strcpy(exprT,"Exact_roots");
        i=save_T(word[4]);

        // 3.4.2017 Coefficient matrix -> Cn.MAT
        C=malloc(m*m*sizeof(double));
        for (i=0; i<m*m; ++i) C[i]=0.0;
        for (i=0; i<m; ++i)
            {
          for (j=0; j<m; ++j)
            C[i+m*j]=T[i]*sgn(cos(pi*T[i+m]*(2*j+odd)/(double)(2*i+1)));
            }
        sprintf(sbuf,"Coefficient_matrix_C%d_for_n=%d",n,n);
        sprintf(name,"C%d",n);

        lab=malloc(8*m);
        numlab(lab,m,8);
        matrix_save(name,C,m,m,lab,lab,8,8,-1,sbuf,0,0);

        external_mat_end(argv1);
        return(1);
        }

iprint(i)
int i;
    {
    sprintf(sbuf," %d",i); sur_print(sbuf); // 13.4.2017
    return(1);
    }

amod(n,i)
int n,i;
    {
    int mod;
    mod=n%i;
    if (mod>(int)(i/2)) mod=i-mod;
    return(mod);
    }

double sgn(x)
double x;
    {
    if (x>0.0) return(1.0);
    if (x<0.0) return(-1.0);
    return (0.0);
    }

double cot(double x)
{ return (1/tan(x)); }


/* _ctest.c 19.5.2017/SM (19.5.2017)
*/

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <survo.h>
#include <survoext.h>

#define MMAX 1001
#define PI 3.141592653589793
unsigned int n,m;
char bin[MMAX];
double  ee[MMAX],rr[MMAX];
int ns[MMAX];
__int64 fs[MMAX];
FILE *temp;

void main(argc,argv)
int argc; char *argv[];
    {
    int i;
    __int64 j;
    double a,accuracy;
    int k,nn;

    if (argc==1) return;
    s_init(argv[1]);
    if (g<2)
        {
        sur_print("\nCTEST n,accuracy                                      ");
        sur_print("\nscans all 2^[(n-1)/2-1] positive linear combinations  ");
        sur_print("\nof edge and chord engths of a n-sided regular polygon ");
        sur_print("\nwith coefficients +1 or -1                            ");
        sur_print("\nand checks that only one of them is equal to          ");
        sur_print("\ncot((2*i+1)*pi/(2*n)), i=1,2,...,(n-1)/2.             ");
        sur_print("\nfor a prime number n.                                 ");
        sur_print("\nExample: CTEST 79,0.0000000000001                     ");
        WAIT; return;
        }
    n=atoi(word[1]);
    m=(n-1)/2;
    accuracy=atof(word[2]);
    for (i=0; i<m; ++i)
       {
       ee[i]=2*sin(((n+1)/2-i-1)*PI/(double)n);
       rr[i]=1.0/tan((2*i+1)*PI/(double)(2*n));
       printf("\nE%d=%g R%d=%g",i+1,ee[i],i+1,rr[i]);
       }

    for (i=0; i<m; ++i) bin[i]='0'; bin[m]=EOS;
    j=(__int64)0;
    nn=0; k=0;

    while (1)
        {
        a=0.0;
        for (i=0; i<m; ++i)
            if (bin[i]=='0') a-=ee[i]; else a+=ee[i];
        ++j;
        ++k; if (k>9999999) { printf("M"); k=0; }
        if (a>0.0) // check positive combinations
            for (i=0; i<m; ++i)
                {
                if (fabs(a-rr[i])<accuracy)
                    { ++ns[i]; fs[i]=j; ++nn; printf("\nX%d\n",nn); }
                }

        // next linear combination: 0 corresponds to coefficient -1
        i=m-1;
        while (bin[i]=='1' && i>=0) --i;
        if (i<0) break;
        bin[i]='1';
        ++i; for (; i<m; ++i) bin[i]='0';
        }

    for (i=0; i<m; ++i)
       if (ns[i]!=1) break;
    if (i<m)
        {
        sur_print("\nNot unique solution for n=%d\n",n);
        for(i=0; i<m; ++i) printf("%d ",ns[i]); getch(); getch(); getch();
        }
    else
        {
        sprintf(sbuf,"CTEST%d.TXT",n);
        temp=fopen(sbuf,"wt");
        fprintf(temp,"%s\n",sbuf);
        for (i=0; i<m; ++i) fprintf(temp,"%I64u\n",fs[i]);
        fprintf(temp,"\n#comb=%I64u",j);
        fclose(temp);
        }
    return;
    }